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Scalar field theory on a lattice falling freely into a 1+1 di- 
mensional black hole is studied using both WKB and numer- 
ical approaches. The outgoing modes are shown to arise from 
incoming modes by a process analogous to a Bloch oscillation, 
with an admixture of negative frequency modes corresponding 
to the Hawking radiation. Numerical calculations show that 
the Hawking effect is reproduced to within 0.5% on a lattice 
whose proper spacing where the wavepacket turns around at 
the horizon is ~ 0.08 in units where the surface gravity is 1. 



I. INTRODUCTION 

If nothing can emerge from a black hole then where 
does the Hawking radiation come from? Relativistic field 
theory tells us that it comes from a trans-Planckian reser- 
voir of outgoing modes at the horizon 10 . The existence 
of such a trans-Planckian reservoir is doubtful on general 
physical grounds. It is the cause of the divergences in 
quantum field theory which are expected to be removed 
in a quantum theory of gravity, and in particular it would 
seem to imply an infinite black hole entropy unless can- 
celed by an infinite negative "bare entropy" . 

String theory has produced an account of the Hawk- 
ing radiation and black hole entropy which requires no 
trans-Planckian reservoir, however in that picture instead 
of a black hole one has a stringy object in flat spacetime 
that arises from a black hole as the string coupling con- 
stant tends to zero [|| . Thus gravitational redshift plays 
no role in the string calculations and nothing is learned 
about the origin of the outgoing modes in a black hole 
spacetimej^] If we could follow the theory as the string 
coupling becomes strong we would presumably learn how 
it is that string theory produces the outgoing black hole 
modes. 

Absent a quantum gravity description of the Hawking 
process, Unruh invented a fluid analog of a black hole ||] 
in which an inhomogeneous flow exceeds the long wave- 
length speed of sound creating a sonic horizon (see also 



) . Quantizing the sound field Unruh argued that the 
horizon would radiate thermal phonons, in analogy with 
the Hawking effect, at a temperature given by h/2ir times 
the gradient of the flow velocity at the horizon. Since 
the sound field has an atomic cutoff, and the physics of 
the fluid is completely understood in principle, the fluid 
model may provide a source of insight into the Hawking 
process. 

The atomic cutoff produces dispersion of sound waves 
leading to subsonic propagation at high frequencies. It 
turns out that this dispersion is all that is needed to ob- 
viate the need for a trans-Bohrian reservoir. In the lin- 
earized theory, an outgoing short wavelength mode can 
be dragged in by the fluid flow and redshifted enough so 
that its increased group velocity overcomes the flow and 
it escapes to infinity. Investigations of a number of linear 
field theory models ]r^-pd| and lattice models p| incorpo- 
rating such high frequency dispersion have demonstrated 
this mechanism of "mode conversion" and have shown 
that the Hawking radiation is recovered]^ The dispersion 
in these models is not locally Lorentz-invariant since a 
frame is picked out in which to specify which frequencies 
are "high" . In the fluid model this is the local rest frame 
of the fluid. A real black hole also defines a preferred 
frame but it does so in a non-local way. Perhaps quan- 
tum gravity produces a dispersion effect related to this 
non-local notion of a preferred frame, or perhaps micro- 
physics is just not locally Lorentz invariant. In any case, 
it seems worthwhile to understand the mechanism of the 
outgoing modes and Hawking radiation in these models 
for the hints it may provide about a correct quantum 
gravity account. 

In PI a falling lattice model was introduced to ad- 
dress the "stationarity puzzle" : how can a low fre- 
quency outgoing mode arise from a high frequency ingo- 
ing mode when it propagates only in a stationary back- 
ground spacetime and hence has a conserved Killing fre- 
quency? In the continuum based dispersive models there 
is no satisfactory resolution of this puzzle. As the outgo- 
ing modes are traced backwards in time their incoming 
progenitors are squeezed into the short distance regime 
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The fact that redshifting plays no role in producing the 
stringy Hawking radiation is presumably tied to the fact that 
the string calculations apply in the "near-extremal" limit that 
the surface gravity tends to zero. 



2 Dispersion can also produce superluminal propagation, 
which allows outgoing modes to emerge from behind the hori- 
zon. The superluminal case does not seem very healthy from a 
fundamental point of view, though it is surely relevant in some 



condensed matter analogues of black hole horizons |13|-|15 
is capable of producing the usual Hawking spectrum fit 



and 
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of the model which is not physically sensible. 

A lattice provides a simple model for imposing a phys- 
ically sensible short distance cutoff, ft might seem at 
first that the most natural choice would be to preserve 
the time translation symmetry of the spacetime with a 
static lattice whose points follow accelerated worldlines. 
On a static lattice, however, the Killing frequency is con- 
served, so outgoing modes arise from ingoing modes with 
the same frequency, and there is no Hawking radiation. 
The in- vacuum therefore evolves to a singular state at the 
horizon. If the lattice points are instead freely falling, a 
discrete remnant of time translation symmetry can still 
be preserved. On such a lattice there is still a frequency 
conservation law, and the unsatisfactory resolution of the 
stationarity puzzle is that the lattice spacing in this case 
goes to zero at infinity (if the lattice points are asymp- 
totically at rest) fl2]| . 

If we insist that the lattice spacing asymptotically ap- 
proaches a fixed constant and that the lattice points are 
at rest at infinity, we find a fascinating resolution of the 
stationarity puzzle [jl2|: the lattice cannot have even a 
discrete time translation symmetry, since there is a grad- 
ual spreading of the lattice points as they fall toward the 
horizon. The timescale of this spreading is of order 1/k 
where n is the surface gravity. This time dependence of 
the lattice is invisible to long wavelength modes which 
sense only the stationary background metric of the black 
hole, but it is quite apparent to modes with wavelengths 
of order the lattice spacing. On such a lattice the long 
wavelength outgoing modes come from short wavelength 
ingoing modes which start out with a high (lattice scale) 
Killing frequency which is exponentially redshifted as the 
wavepacket propagates to the horizon and turns around. 

This behavior of wavepackets in the (time-dependent) 
falling lattice model was determined in [O with the help 
of the WKB approximation in which the frequency is 
used as a Hamiltonian to solve for the wavepacket trajec- 
tory. It was also argued that, since the time dependence 
of the lattice is adiabatic for those modes with wave- 
length short enough to know they are on a lattice, these 
modes will remain unexcited on their way to the horizon 
and the ground state condition for the Hawking effect 
will be met. However, since the wavelength is of order 
the lattice spacing and the wavepacket can vary wildly 
on the lattice near the horizon, it is not obvious that the 
WKB approximation is reliable nor that the adiabatic 
argument really applies. 

The primary motivation for the work reported here 
was to check the WKB and adiabatic assumptions by 
carrying out the exact calculation numerically using the 
lattice wave equation. We found that indeed the field be- 
haves this way, and the Hawking radiation is recovered 
to within half a percent for a lattice spacing S = 0.002/k 
which corresponds to a proper spacing ~ 0.08/k where 
the wavepacket turns around at the horizon. We also 
studied the deviations from the Hawking effect that arise 
as k8 is increased, and our simulations reveal an inter- 
esting picture of how the wavepackets turn around at the 



horizon. 

The precise choice of worldlines for the lattice points 
should not be important to the leading order Hawking 
effect as long as the ground state evolves adiabatically in 
the lattice theory. Our numerical results are consistent 
with this expectation, in that the particular worldlines 
depend on the lattice spacing yet the results converge to 
the continuum as the lattice spacing is decreased. More- 
over it was shown recently JTcj ] , in a continuum model 
with high frequency dispersion, that when the preferred 
frame is changed from the free-fall frame of the black 
hole to that of a conformally related metric there is no 
leading order change in the Hawking radiation unless the 
acceleration of the preferred frame is very large. Pre- 
sumably all that matters is that the preferred world lines 
flow smoothly across the horizon. From a fundamental 
point of view, if there really is a cutoff in some preferred 
frame, it is plausible that this frame would coincide with 
the cosmic rest frame and would fall across the event 
horizons of any black holes that form by collapse. 

The remainder of this paper is organized as follows. In 
Section |n] the falling lattice model is set up, in Section 
III some results of the WKB method are shown, illustrat- 
ing the frequency shift and the mode conversion at the 
horizon, and the role of the negative frequency branch is 
explained. In Section |^ the results of the full numer- 
ical evolution are presented, and Section ^ contains a 
discussion or the results and directions for further work. 
The Appendix describes the finite differencing of the time 
derivates used in the numerical calculations. 

We use units with K = c = k = 1, where k is the 
surface gravity. 



II. FALLING LATTICE MODEL 

In this section we set up the field theory on a lattice 
falling into a black hole in two spacetime dimensions. 

A. Black hole spacetime 

We assume the spacetime metric is static, so [[l8| co- 
ordinates can be chosen (at least locally) such that the 
line element takes the form 

ds 2 = dt 2 ~ (dx - v(x)dt) 2 . (2.1) 

In these coordinates the Killing vector is given by 

X = d t , (2.2) 

with squared norm 

X 2 = l~v 2 {x). (2.3) 
For v{x) we choose 



v(x) 



cosh(/3a;) 



(2.4) 
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where 



and (3 are positive constants. As \x\ 



+00, v(x) vanishes, so the line element becomes that of 
Minkowski spacetime. Provided v max > 1 there is an er- 
goregion centered on x = in which the Killing vector 
is spacelike, and the boundaries of this region are black 
and white hole horizons at 



x H = ±(i 1 cosh l {v max ). 
The surface gravity of the horizons is given by 



k = v'(xh) = ft vl V; 



-2 
max- 



(2.5) 



(2.6) 



The surface gravity sets the length scale for our space- 
time, and we will keep it fixed. Hence it will be useful to 
employ units for which k = 1. 

We now introduce a freely falling "Gaussian" coordi- 
nate z which will be discretized to define a falling lattice. 
The worldlines with dx — v(x)dt are geodesies of the 
metric (2.1) with proper time t, at rest as \x\ — ► 00, 
and they are orthogonal to the surfaces of constant 
t. On these curves the quantity t — Jdx/v is con- 
stant, so if z is constant on these curves we must have 
W(z(x, t j) — t — J dx/v for some function W. We choose 
this function so that z = x at t = 0, which implies that 
W(z) = — J z dx/v = sinh(/3z)//3i; max , hence 



sinh(/3z) = sinh(/?a;) + f3v max t. 



(2.7) 



In terms of z, the line element (2^1) takes the form 

ds 2 = dt 2 - a 2 (z,t)dz 2 , (2. 



with 



a(z, t) 



cosh(/3z) 



yjl + (sinh(/3z) - (3v max t¥ 



(2.9) 



We call a(z,t) the "scale factor" , even though it depends 
on z as well as t. It can also be expressed as 



a(z, t) — 



v(x{z, t)) 
v(z) 



(2.10) 



In (t, z) coordinates the Killing vector (2.2) is given by 

(2.11) 



X = d t - v{z)d z =d t ; — —d z 



a(z, t) 

The scale factor is unity at t = 0, and outside the black 
hole it grows towards the future. It is important to know 
how rapidly the scale factor is changing in time, since in 
the lattice theory that time dependence is not merely a 
coordinate effect and can therefore excite the quantum 
field. At the horizon the scale factor takes the value 

a(z H , t)^^l + 2nt + kH 2 /{1 - v m L), (2.12) 

so a(zH,t) Hd for nt > 1. Thus the fractional rate of 
change dta/a is of order n. In general, we have 



dta/a 



f3v max smh((3x) 
cosh 2 (/3a;) 



(2.13) 



whose maximum outside the horizon (cosh(f3x) > w max ) 
is k (at the horizon) if u max > V% and CTmax/V^max — 1 
if Umax < \/2. Thus \dta/a\ < O(k) everywhere outside 
the horizon as long as u max > 1 is not too close to unity. 



B. Lattice 

The lattice is defined by discretizing the z coordinate 
with spacing Az = S. That is, we have only the discrete 
values z m = mS. We will consider values of k8 which 
range from 0.002 to 0.128. To give a feeling for the lattice 
we plot in Fig. [j] the (t, x) coordinates of lattice point 
worldlines at intervals of 50 lattice points. 
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FIG. 1. Every 50th lattice point, t vs. x, for the case 
5 — 0.004 and i> max = v2- In the black region of the plot 
the points are too close to resolve. The horizon is at x ~ 0.6. 

Since dx — a{z 1 t)dz at constant t, the separation Aa; 
of the lattice points in Fig. [l] gives a direct representation 
of the scale factor. The prop er lat tice spacing is given by 
a(z,t)5 which, according to ( |2.12 ), grows approximately 
linearly with time at the horizon. 

Note that there is a region where the lattice spacing 
becomes very small. This happens because we chose the 
spacing to be uniform at t = 0. If the points did not 
bunch up outside the black hole before t = 0, they could 
not wind up equally spaced at t = Q, since they are all 
on unit energy free-fall trajectories. (The line element is 
symmetric under (t, x) — » (— t, — x), so the points bunch 
up also outside the white hole horizon after t — 0.) 

To get a better idea of how sparse the lattice is near 
the horizon at the turning point of a typical wavepacket, 
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we plot in Fig. |] the velocity function v (2.4) for every 
lattice point near the horizon at time t = 27, for the 
case u m ax = V% and kS — 0.004. In this case there 
are only 6 points in the region between v(x) = 0.5 and 
the horizon, and the proper spacing at the horizon is 
~ y/2nSt ~ 0.15/ 
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FIG. 2. v(x) for every lattice point at t — 27, for the case 
Umax = s/2 and S = 0.004. 



C. Scalar field 

We consider a massless, minimally coupled scalar field 
on the lattice. In the continuum the action would be 



Scant — 



dtdz (a(z, t)(d t ipY 



a(z, t) 



(c^) 2 . (2.14) 



where we substituted the metric ( |2.S| ) in the second line. 
On the lattice we adopt the same action except that the 
partial derivative d z ip is replaced by the finite difference 
Dip m := (cpm+i — (fm)/S and a m (t) is replaced by its 
average over the relevant lattice sites: 

2 J \ a m+ i(t) + a„ 



(2.15) 



Varying the action ( 2.1 5| ) gives the equation of motion 

for tp m (t), 

d t (a m (t)d t cp m (t)) - D ( 2D f" l - l{t) ^ ) = 0. (2.16) 



(t)+a m (t) 

There is a conserved (not positive definite) inner product 
between complex solutions: 

(Vw) = i^a m {t){il>* m {t)dt(p m {t) - (p m (t)d t ^{t)). 

m 

(2.17) 



In two dimensions the continuum action (2.14) is con- 
formally invariant, and all metrics are conformally flat, 



hence there is no wave scattering. The discretization 
breaks this conformal invariance, however there is still 
essentially no scattering on the lattice providing the lat- 
tice spacing is much smaller than the radius of curvature. 
The reason is that modes with wavelength long compared 
to the lattice spacing cannot tell they are not in the con- 
tinuum, while modes with wavelength comparable to the 
lattice spacing cannot tell that the spacetime is curved. 
Thus, when we compute the Hawking occupation num- 
bers, there is no "greybody factor" . 



D. Quantization and the Hawking effect 

The field is quantized as a collection of self-adjoi nt op - 
erators m (t) satisfying the equation of motion ( 2.16| ) 
and canonical commutation relations. For each complex 
solution to the equations of motion we defin e the oper- 
ator a(f) — (f)0) using the inner product ( 2.17 ). The 
commutation relations are equivalent to the relations 



W),ai(g)] = (J,g) 



(2.18) 



for all / and g. 

For a positive norm solution p, a(p) acts like a lowering 
operator, while for a negative norm solution n, the con- 
jugate n* has positive norm, so a(n) = —a^(n*) acts like 
a raising operator. The Hilbert space of ingoing modes 
is just the Fock space built from solutions with positive 
frequency with respect to t (hence positive norm), and 
similarly for the outgoing Hilbert space. We assume the 
boundary condition that the ingoing positive frequency 
field modes are all in their ground state, 



a(p in )\0) = 0. 



(2.19) 



The occupation number (0\aJ (p)a(p)\0) of an outgo- 
ing normalized positive frequency wavepacket p in the 
in- vacuum |0) is is just minus the norm of the negative 
frequency part of the ingoing wavepacket that gives rise 
to p. To see this, suppose the ingoing wavep acket q+ + q- 
evolves to p according to the field equation (2.16), where 
q+ has positive frequency and q- has negative frequency. 
Then a(p) = a(q+) + a(q_) = a(q+) — a^gl), hence 



(0|ot(p)o(p)|0) =-(?_, 



(2.20) 



where we used the in-vacuum condition (2.19) and the 
commutation relations (2.15). If p is not normalized then 
to obtain the occupation number we must divide by (p, p). 
Hawking's result |13] in standard field theory yields, for a 
wavepacket p = fdw c u exp(~ioj(k)t + ikz), the thermal 
occupation number 



-^Hawking = / Cfcj 



J/Tn - 1 ■ 



(2.21) 



where Th = k/2tt is the Hawking temperature. 
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III. WKB ANALYSIS 

The motion of wavepackets in the falling lattice model 
can be studied using the WKB approximation [Q . This 
amounts to using Hamilton's equations for the position 
z and momentum k with a Hamiltonian determined by 
the dispersion relation 



H(z,k,t) =uj = ± 



2/J 
a(z, t) 



sin(fe<ty2) 



(3.1) 



which is obtained by inserting a mode of the form 
exp(— iu>t + ikz) into the held equation and keeping only 
the terms with the highest derivatives or finite differences 
of the field . A plot of this dispersion relation is shown 
in Fig. H for a(z,t) = 1. Wavevectors differing by 2tt/S 
are equivalent, so only the range (—ir/S,ir/5) (the "Bril- 
louin zone") is plotted. 
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FIG. 3. Dispersion relation uj8 = ±2 sin(fc<5/2) plotted 
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The frequency u) is not conserved, even when kS -C 1, 
simply because it is not the Killing frequency. Using 
(2.11) we see that in the continuum the Killing frequency 



would be given by 



L0 K 



v(x(z,t)) 
OJ H ; — - — k, 



(3.2) 



a(z, t) 

with which the dispersion relation can be expressed as 

a(z, t)u K ~ v{x)k = (2/5) sin(fc<5/2). (3.3) 

Note however that this relation is only valid when kS 
is not too large, otherwise the fact that the Killing fre- 
quency involves a partial derivative in a direction that is 
not along a lattice point worldline renders the relation 
meaningless. 

It was shown in JlJ] that the WKB trajectory of an in- 
going positive frequency wavepacket bounces off the hori- 
zon and comes back out provided the ingoing wavevec- 
tor is greater in magnitude than some critical value \k c \. 
(This critical value depends on the time and place from 
which the packet is launched.) The wavevector evolves 
from negative values to positive values by decreasing un- 
til it is less than —tt/S, at which point the group velocity 



(in the free fall frame) changes sign and the wavevector is 
equivalent to +tt/S. From there it continues to decrease, 
winding up small and positive. This reversal of group 
velocity by a monotonic change of the wavevector on a 
lattice is analogous to the Bloch oscillation of an electron 
in a crystal acted on by a uniform electric field. 

Thus the interval (— n/S, — \k c \) is mapped onto the in- 
terval (0,7r/5). The map is onto since every outgoing 
wavevector arises from some ingoing wavevector. This 
dynamical stretching of the fe-space interval of length 
tt/S — | fe c | to one of length tt/S is not problematic in a 
particle phase flow, but if we think about the wavepack- 
ets that are being approximated by this WKB analysis it 
is disturbing since there really are more outgoing modes 
than ingoing modes which produced them. It seems the 
wave evolution is somehow not conserving the number of 
wave degrees of freedom. 

The resolution of this puzzle is the essence of the Hawk- 
ing effect. The wave evolution mixes positive and nega- 
tive frequency modes, so a purely positive frequency out- 
going wavepacket arises from a mixture of positive and 
negative frequency ingoing wavepackets. WKB evolu- 
tion breaks down at the turing point, hence it misses this 
mode mixing. 

An example of a positive frequency WKB trajectory is 
shown in Fig. ^. The graph shows several different quan- 
tities as functions of t (in units with k = 1): the static 
position coordinate of the wavepacket x an d th e horizon 
xh (solid lines), the free-fall frequency lj ([O]) (dashed 
line), the Killing frequency uik ( 3~2] ) (dot-dashed line), 
and 50sin(fe<5/2) (sparsely-dashed line). The positions 
are scaled up by a factor of 10 so they can be seen on 
the same graph. In this example the lattice spacing is 
8 = 0.004, w max = V2, and the trajectory has z(34) = 6 
and fe(34) = 1. 
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FIG. 4. WKB evolution. 



Notice first the behavior of sin(fc<5/2). Backwards in 
time this starts out very small, rises to the maximum, 
and comes back down but not to zero. The group veloc- 
ity in the frame falling with the lattice vanishes at the 
top, when kd/2 = tt/2, corresponding to a wavelength 
equal to two lattice units. This rise and fall happens 
via a monotonic change of k, and is the Bloch oscillation 
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mentioned above. 

The wavepacket that ends up with u = ll>k = 1 at 
late times starts out at early times with u) — ujk ~ 43. 
The maximum possible lattice frequency is 2/5 = 500 
in this case, so the ingoing frequency is certainly "high" . 
Whereas the free-fall frequency continues to redshift (for- 
ward in time) as the wavepacket climbs away from the 
horizon, the redshifting of the Killing frequency is essen- 
tially complete by the time the wavepacket reaches the 
turning point. Thus, past the turning point, each Killing 
frequency component propagates autonomously, more or 
less as it would in the continuum. As described in the 
Introduction, this, together with the fact that the evolu- 
tion of u> at the lattice scale is adiabatic, is the reason 
why the usual Hawking effect is reproduced. 



IV. LATTICE WAVE EQUATION ANALYSIS 

In this section we discuss the results of our numerical 
calculations. The basic method is to take an outgoing 
positive frequency wavepacket and evolve it backwards 
in time. After it has "bounced" from the horizon and 
propagated back into the flat region of the spacetime we 
decompose it into positive and negative frequency parts. 
Since it is a purely left moving wavepacket at this stage 
the negative frequency part is just the part com posed 
of positive wavevectors. The occupation number ( 2.20| ) 
of the original outgoing packet is then given by minus 
the norm of this negative frequency part divided by the 
total norm. We typ ically computed the total norm in 
position space using (2.17). Since the negative frequency 
part if(~^ was projected out in the wavevector space we 



used an alternate expression for its norm, [tp^ ' , <p( )) 



-2E fe >oK fc )P M (fc)l : 

given by (£ 



where the frequency u>(k) is 



A. Numerical issues 



Our system of equations ( 2 . 1 6| ) is a coupled set of or- 
dinary differential equations, which is to be solved in the 
limit that time is continuous but the spatial lattice is 
fixed. In order to satisfy the Courant stability criterion 
the time step At must satisfy At < a(z, t)S. (In this 
1+1 dimensional setting the Courant condition coincides 
with the condition that causality be respected by the dif- 
ferencing scheme.) In the asymptotic region a(z, t) = 1 
so the stability condition there is At < 5. For a fixed 
At/5 the general condition would be violated wherever 
a(z, t) < At/5. Since this only happens deep behind the 
horizon (for t > 0) we just modified a(z, t) inside the 
black hole so that it is never less than At/ 5. We used a 
time discretization scheme (written out in the Appendix) 
with errors of order (At) 2 , and found that At = QA5 was 
adequately small to obtain very good accuracy. 



B. Wavepackets and parameters 

If the Killing frequency were conserved as it is in 
the dispersive continuum models the entire Hawking 
spectrum could be probed just by propagating one 
wavepacket and decomposing it into is Fourier compo- 
nents, each of which would pick up a negative frequency 
part appropriate for the corresponding frequency. This 
is how Unruh checked the spectrum in his model [p|. 
Since the Killing frequency is not conserved on our lattice 
we cannot avail ourselves of this option.^] Thus instead 
we compute the occupation numbers for a sequence of 
wavepackets Lp s (k) with different profiles: 



<p B {k) = k exp[-(fc - 0.05s) 2 ], 



s = l, 



,9 (4.1) 



(in units with k = 1, as usual). The form of these 
wavepackets was dictated by the need to have them well 
enough localized to be contained in the flat region of 
the spacetime (so we could construct the correspond- 
ing positive frequency initial data) while at the same 
time containing enough power at low frequencies of or- 
der k to have a measurable Hawking occupation num- 
ber. To shift the wavepacket t o th e desired starting po- 
sition the Fourier transform (4.1) was multiplied by a 



phase factor exp(— i/czimtiai)- The positive frequency ini- 
tial data was constructed by evolving each Fourier com- 
ponent one time step At by multiplying with the phase 
factor exp(— iu>(k) At) and Fourier transforming the re- 
sult back into position space. 

The time dependence of the lattice (illustrated in Fig. 
|l|) means that the results will not be independent of 
the "launch time" of the wavepacket from a given lo- 
cation. If we start our backwards evolution too far in 
the future then when the wavepacket nears the horizon 
the lattice points will be pathologically sparse. On the 
other hand if we launch backwards at too early a time 
then the wavepacket can become partially trapped in the 
region where the lattice spacing gets very small. This 
happens because a sufficiently high frequency wave is 
turned around when it tries to propagate into a region 
of larger lattice spacing. We encountered this "chan- 
neling" phenomenon in our simulations, and avoided it 
just by launching later. The runs reported here were all 
done with wavepackets that reached the horizon around 
t = 27, and were launched from far enough away to be 
contained in the flat region of the spacetime. Typically, 



We studied the possibility that the WKB evolution could 
be used to establish a mapping between in and out frequen- 
cies, thus allowing us to check the entire spectrum of particle 
creation just by propagating a single wavepacket. This turns 
out to be a flawed idea however since the modes are mixed 
by the evolution and, moreover, such a map misses the role 
of the negative frequency piece. 
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the launch time was around t = 56 when the wavepacket 
was centered at x = 28. 

All of the calc ulat ions reported her e were done using 
the line clement (2.1), with u max (2.4) ranging from \/2 
to 30. 



A sample of the comparisons for kS = 0.002 is listed 
in Table 3. The typical relative difference in this case is 
of order half a percent. 



3. Lattice dependence 



C. Results 



1. Generalities 



The behavior of a typical wavepacket throughout the 
process of bouncing off the horizon is illustrated in Fig. 

The real part of the wavepacket is plotted vs. the 
static coordinate x at several different times. Following 
backwards in time, the wave starts to squeeze up against 
the horizon and then a trailing dip freezes and develops 
oscillations that grow until they balloon out, forming into 
a compact high frequency wavepacket that propagates 
neatly away from the horizon backwards in time. 

This ingoing wavepacket contains both positive and 
negative frequency components, in just the right combi- 
nation to produce only an outgoing wave when sent in 
towards the black hole, with no wave propagating across 
the horizon. To illustrate this we have taken one of these 
ingoing wavepackets, decomposed it into its positive and 
negative frequency parts, and sent them separately back 
toward the horizon forward in time (Fig. |7]). The nega- 
tive frequency part (which is the smaller of the two since 
the surface gravity is fairly low compared to the typical 
frequencies in the outgoing wavepacket), mostly crosses 
the horizon, with just a bit bouncing back out. The pos- 
itive frequency part mostly bounces, with just a little bit 
going into the black hole. 



2. Spectrum 

T he " observed" occupation numbers for the wavepack- 
ets (4.1) all a gree well with the integrated thermal pre- 
dictions ( 2.21 ) of the Hawking effect provided the lattice 
spacing is not too large. 



TABLE I. Comparison of the black hole radiatio n on the 
lattice to the thermal Hawking occupation numbers (2.21) for 
the wavepackets (1.1) for the case i> max — V2 and k8 = 0.002. 



s 


Thermal 


Lattice 


Rel. Diff. 


1 


0.01563 


0.01557 


0.0038 


2 


0.01409 


0.01404 


0.0035 


3 


0.01266 


0.01262 


0.0034 


4 


0.01135 


0.01130 


0.0041 


5 


0.01014 


0.01009 


0.0045 


6 


0.00903 


0.00899 


0.0039 


7 


0.00801 


0.00797 


0.0051 


8 


0.00709 


0.00705 


0.0052 


9 


0.00625 


0.00622 


0.0046 



To study the dependence on the lattice we ran the 
s = 1 wavepacket backward at several different valu es o f 
the lattice spacing and for several values of v max (2.4) 
from v2 to 30. The change of u max produces only a 
rather small change in the shape of the function v(x) 
outside the horizon, ho wev er it produces a shift of the 
position of the horizon (2.5). For large w max the horizon 
is given approximately by xh — hr(2t> max ). The horizon 
position ranges from about 0.6 for t> max = \[2 to about 4 
for w max = 30. Thus a wavepacket launched backward in 
time from a given position at a given time will reach the 
horizon later in time for larger w max , so that the proper 
lattice spacing at the horizon will be larger, though not 
by a very large amount for the wavepackets considered 
here. 

In Fig. U we plot the occupation numbers for lattice 
spacings 5 rang ing from 0.002 to 0.128, and for several 
values of v. 



«W = V2, 2, 2.5, 3, 4, 5, 10, 20, 30. 
In all cases the late time wavepacket was launched from 
the same position and time. All wavepackets had the 
same envelope for their Fourier transform, however for 
different lattice spacings the wavepacket is necessarily 



different. 
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FIG. 5. Occupation numbers for s = 1 wavepacket for var- 
ious values of 5 and w max . Higher occupation numbers corre- 
spond to larger values of u max . 

For all values of u max , the occupation number con- 
verges to the Hawking prediction as the lattice spacing 5 
decreases. 

The pattern of deviations from the thermal prediction 
appears to depend quite strongly on the value of f max , a 
result which we do not understand at present. The best 
we can do is to list various effects which would contribute 
to the deviations: (i) the WKB turning point is mov- 
ing away from the horizon as S grows, so one might ex- 
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pect that the effective surface gravity v'(x t . p .) and hence 
the effective Hawking temperature felt by the wavepacket 
would decrease; (ii) at larger u max the proper lattice spac- 
ing at the horizon is larger at the time this particular 
wavepacket reaches the horizon; (iii) the lattice points 
are rather sparse near the horizon for the larger values 
of S, and their precise positions will be significantly dif- 
ferent for different values of v max . Effect (i) would lead 
to a decrease in the occupation number, while one might 
expect that effect (ii) would lead to a relative increase 
due to the extra particle creation associated with time 
dependence of the lattice. Perhaps what we are seeing is 
just a delicate balance between these two effects. 

V. DISCUSSION 

We have shown that linear field theory on a falling lat- 
tice reproduces the continuum Hawking effect with high 
fidelity, thus verifying earlier expections |Q. For the 
smallest lattice spacing we studied, S = 0.002/k, the de- 
viations from the thermal particle creation in wavepack- 
ets composed of frequencies to < O(k) were of order 
half a percent. This is remarkable, particularly since the 
proper lattice spacing encountered by the wavepacket at 
the horizon was ~ 0.08/k. When kS is pushed to higher 
values eventually we see significant deviations from the 
Hawking prediction. 

There are several directions in which this work could 
be developed. One idea is to include self-interactions for 
the quantum field theory. There are general arguments 
pl| , p2| that the thermal Hawking effect holds also for in- 
teracting fields. There are also perturbative calculations 
|p3| and calculations exploiting exact conformal invari- 
ance of certain 1+1 dimensional field theories Also 
the "hadronization" of Hawking radiation from primor- 
dial black holes has been analyzed by matching to high 
energy collision data [^5| , but the physics of how the ther- 
mal Hawking radiation is "dressed" by interactions as it 
redshifts away from the horizon has not been studied. In 
a 1+1 dimensional model it should not be prohibitively 
difficult to numerically study Hawking radiation of inter- 
acting fields, now that we see it can be done on a lattice. 

The falling lattice model has provided a satisfactory 
mechanism — the Bloch oscillation — for how to get an 
outgoing mode from an ingoing mode in a stationary 
background, but there is a serious flaw in the picture: 
the lattice is constantly expanding. In the fluid model by 
contrast the lattice of atoms maintains a uniform average 
density. In a fundamental theory we might also expect 
the scale of graininess of spacetime to remain fixed at, 
say, the Planck scale or the string scale (since presum- 
ably the graininess would define this scale) rather than 
expand. Can the falling lattice model be improved to 
share this feature? 

A fluid maintains uniform density in an inhomogeneous 
flow by compressing in some directions and expanding 



in others, which requires at least two dimensions to be 
possible. At the atomic level such a volume-preserving 
flow involves erratic motions of individual atoms. One 
possible improvement of the lattice model is to make a 
lattice that mimics this sort of volume preserving flow. It 
is not clear whether the motions of the lattice points can 
be slow enough to be adiabatic on the timescale of the 
high frequency lattice modes. If they cannot, then the 
time dependence of this erratic lattice background will 
excite the quantum vacuum. 

In a fluid the lattice is a part of the system, not just 
a fixed background. The average flow could be adiabatic 
for the fully coupled ground state of the system but not 
for the field theory of the perturbations on the "back- 
ground lattice" . Similar comments apply in quantum 
gravity: surely if in-out mode conversion is at play the 
incoming high frequency modes are strongly coupled to 
the quantum gravitational vacuum. Ideally, therefore, we 
should try to find a model in which the background is not 
decoupled from the perturbations. 

As a first step in this direction one could study a one 
dimensional model in which the lattice points are non- 
relativistic point masses, coupled to each other by near- 
est neighbor interactions, and "falling" or propagating in 
a background potential (with or without periodic bound- 
ary conditions). The perturbations of such a lattice are 
the phonon field, and the back reaction to the Hawking 
radiation is included (although the background potential 
is fixed). In a model like this one could presumably follow 
in detail the nonlinear origin of the outgoing modes and 
the transfer of energy from the mean flow to the thermal 
radiation. 
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APPENDIX: FINITE DIFFERENCING SCHEME 



The equation of motion ( 2.16 ) for ip m (t) is 

dt{a m (t)d t f m {t) 



= D 



2£><p ro _ 1 (i) 



(t)+a m (t) / 

We used a time discretization scheme for the left hand side arising from the replacement of / by 
>(t + At/2) - f(t - At/2)] /At = f+ (l/24)/( 3 )(Ai) 2 + 0((Atf), 

which yields 

dt(a m (t)d t <p m {t)) -» [o(n + 1/2) (v»(n + 1) - p( n )) - a (« - V 2 ) ~ - !))] /(Ai) 2 , 



(Al) 



(A2) 



(A3) 



where the spatial index is omitted and the argument n stands for nAt. We could have used this replacement, since 
a{z, t) is a prescribed function that can be calculated at any t, but for some reason we substituted a(n ± 1/2) by the 

average ( a(n) + a(n ± 1)J /2. This yields the approximation 



a(n) + a(n + 1)J Up{n + 1) - p(n)J - fa(rc) + a(n - 1)) f<p(n) - y>(n - 1)J /2(At) 2 
= 9 t (aS t <^)(n) + (^ (4) + ^ (3) + + J« (3 V) H(Ai) 2 + 0((Ai) 4 ). 



(A4) 
(A5) 



With (A4) in place of the left hand side of (Al) we solve for ip(n — 1) in terms of ip(n) and (p(n + 1) to obtain an 
equation for evolving ip backwards in time by one time step given the following two steps. 
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FIG. 6. A typical wavepacket evolution. Here tWx = \/2, k<5 = 0.004, and the wavepacket ( J4. l[ ) has s = 1. The ocsillations 
of the incoming wavepacket are too dense to resolve in the plots. 
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t=6 : = 30 t = 54 




FIG. 7. The positive (top row) and negative (bottom row) frequency parts of the ancestor of an outgoing positive frequency 
wavepacket, sent back towards the black hole. Note the vertical scale is different in the two plots. The wavepackets crossing the 
horizon in the two cases are negatives of each other since they must cancel to produce only the outgoing wavepacket we started 
with. The parts that go through the horizon appear on the right due to periodic boundary conditions. The high frequency 
ocsillations are too dense to resolve in the plots. 
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